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Abstract. Here we introduce the Delaunay Density Es- 
timator Method. Its purpose is rendering a fully volume- 
covering reconstruction of a density field from a set of 
discrete data points sampling this field. Reconstructing 
density or intensity fields from a set of irregularly sam- 
pled data is a recurring key issue in operations on as- 
tronomical data sets, both in an observational context as 
well as in the context of numerical simulations. Our tech- 
nique is based upon the stochastic geometric concept of 
the Delaunay tessellation generated by the point set. We 
shortly describe the method, and illustrate its virtues by 
means of an application to an N-body simulation of cos- 
mic structure formation. The presented technique is a fully 
adaptive method: automatically it probes high density re- 
gions at maximum possible resolution, while low density 
regions are recovered as moderately varying regions devoid 
of the often irritating shot-noise effects. Of equal impor- 
tance is its capability to sharply and undilutedly recover 
anisotropic density features like filaments and walls. The 
prominence of such features at a range of resolution levels 
within a hierarchical clustering scenario as the example of 
the standard CDM scenario is shown to be impressively 
recovered by our scheme. 

Key words: Methods: numerical, statistical, N-body sim- 
ulations - Cosmology: large-scale structure of the Universe 



1. Introduction 

Astronomical observations, physical experiments as well 
as computer simulations often involve discrete data sets 
supposed to represent a fair sample of an underlying 
smooth and continuous field. Conventional methods are 
usually plagued by one or more artefacts. Firstly, they of- 
ten involve estimates at a restricted and discrete set of 
locations - usually defined by a grid - instead of a full 
volume-covering field reconstruction. A problem of a more 
fundamental nature is that the resulting estimates are im- 
plicitly mass-weighted averages, whose comparison with 
often volume-weighted analytical quantities is far from 
trivial. For most practical purposes, the disadvantage of 
almost all conventional methods is their insensitivity and 



infiexibility to the sampling point process. This leads to 
a far from optimal performance in both high density and 
low density regions, which often is dealt with by rather 
artificial and ad hoc means. 

In particular in situations of highly non- uniform distri- 
butions conventional methods tend to conceal various in- 
teresting and relevant aspects present in the data. The cos- 
mic matter distribution exhibits conspicuous features like 
filaments and walls, extended along one or two directions 
while compact in the other(s). In addition, the density 
fields display structure of varying contrasts over a large 
range of scales. Ideally sampled by the data points, appro- 
priate field reconstructions should be set solely and auto- 
matically by the point distribution itself. The commonly 
used methods, involving artificial filtering through for in- 
stance grid size or other smoothing kernels (e.g. Gaussian 
filter) often fail to achieve an optimal result. 

Here we describe and propagate a new fully self- 
adaptive method based on the Delaunay triangulation of 
the given point process. After a short description of the 
fundamentals of our tessellation procedure, we show its 
convincing performance on the result of an N-body simu- 
lation of structure formation, whose particle distribution 
is supposed to reflect the underlying cosmic density field. 
A detailed specification of the method, together with an 
extensive quantitative and statistical evaluation of its per- 
f ormance will be presented in a fort hcoming publication 
( gchaap and van de Weygaert, 200(j ). 



2. The Delaunay Tessellation Field Estimator 

Given a set of field values sampled at a discrete number 
of locations along one dimension we are familiar with var- 
ious prescriptions for reconstructing the field over the full 
spatial domain. The most straightforward way involves 
the partition of space into bins centered on the sampling 
points. The field is then assumed to have the - constant 
- value equal to the one at the sampling point. Evidently, 
this yields a field with unphysical discontinuities at the 
boundaries of the bins. A first-order improvement con- 
cerns the linear interpolation between the sampling points, 
leading to a fully continuous field. 



2 



W.E. Schaap & R. van de Weygaert: Delaunay Reconstruction of Continuous Fields 





/ 1 




Y ..•■'/ 'v,^ /*\" 






"■ k'l 


/ '^'"^^ \ \ s - - / 

/ , . - K Ni / 


- V 4f.r 


_ J •■7 / y >? 

/ \ "vs: " \ 1 

■•■■■■jf- " V---''' / \ 






Fig. 1. A set of 20 points with their 
Voronoi (left frame: solid lines) and De- 
launay (right frame: solid lines) tesse- 
lations. Left frame: the shaded region 
indicates the Voronoi cell correspond- 
ing to the point located just below the 
center. Right frame: the shaded region 
is the "contiguous Voronoi cell" of the 
same point as in the lefthand frame. 



In more than one dimension, the equivalent spatial in- 
tervals of the 1-D bins are well-known in stochastic ge- 
ometry. A point process defines a Voronoi tessellation by 
dividing space into a unique and volume- covering network 
of mutually disjunct convex polyhedral cells, each of which 
comprises that part of multidimensional space closer to 
the defining point than to any of the other (see van de 
Weygaert (1991) and references therein). These Voronoi 
cells (see Fig. 1) are the multidimensional generalization 
of the 1-D bins in which the zeroth-order method approxi- 
mates the field value to be constant. The natural extension 
to a multidimensional linear interpolation interval then 
immediately implies the corresponding Delaunay tessella- 



tion (Delone, 1934). This tessellation (Fig. 1) consists of 
a volume-covering tiling of space into tetrahedra (in 3-D, 
triangles in 2-D, etc.) whose vertices are formed by four 
specific points in the dataset. The four points are uniquely 
selected such that their circumscribing sphere does not 
contain any of the other datapoints. The Voronoi and De- 
launay tessellation are intimately related, being each oth- 
ers dual in that the centre of each Delaunay tetrahedron's 
circumsphere is a vertex of the Voronoi cells of each of 
the four defining points, and conversely each Voronoi cell 
nucleus a Delaunay vertex (see Fig. 1). The "minimum tri- 
angulation" property of the Delaunay tessellation has in 
fact been well-known and abundantly applied in, amongst 
others, surface rendering applications such as geographical 
mapping and various computer imaging algorithms. 

Consider a set of N discrete datapoints in a finite re- 
gion of M-dimensional space. Having at one's disposal 
the field values at each of the (1-t-M) Delaunay vertices 
Xq, Xi, . . . , xm, at each location x in the interior of a De- 
launay M-dimensional tetrahedron the linear interpolation 
field value is defined by 



/(x) = /(xo) + V/(xo)|d,i(x-xo), 



(1) 



in which V/(xo)|]3(,j is the estimated constant field gradi- 
ent within the tetrahedron. Given the (1-l-M) field values 
/(xg), /(xi), . . . , /(xm), the value of the M components of 
V/(xo)|j3(,j can be computed straightforwardly by evalu- 
ating Eqn. (1) for each of the M points Xi, . . . , Xm- This 
multidimensional procedure of linear interpolation was al- 



ready described by Bernardeau & van de Weygaert ( |1996D 
in the context of defining procedures for volume- weighted 
estimates of cosmic velocity fields. While they explicitly 
demonstrated that the zeroth-order Voronoi estimator is 
the asymptotic limit for volume- weighted field reconstruc- 
tions from discretely sampled field values, they showed the 
superior performance of the first-order Delaunay estima- 
tor in reproducing analytical predictions. 

The one factor complicating a trivial and direct imple- 
mentation of above procedure in the case of density (inten- 
sity) field estimates is the fact that the number density of 
data points itself is the measure of the underlying density 
field value. Unlike the case of velocity fields, we there- 
fore cannot start with directly available field estimates at 
each datapoint. Instead, we need to define appropriate es- 
timates from the point set itself. Most suggestive would 
be to base the estimate of the density field at the location 
Xi of each point on the inverse of the volume Vvor,i of its 
Voronoi cell, p{xi) = m/Vvor,i- Note that in this we take 
every datapoint to represent an equal amount of mass m. 
The resulting field estimates are then intended as input 
for the above Delaunay interpolation procedure. However, 
one can demonstrate that integration over the resulting 
density field would yield a different mass than the one 
represented b y the set of sample points (see Schaap & van 
de Weygaert ( ^OOOD for a more specific and detailed discus- 
sion). Instead, mass conservation is naturally guaranteed 
when the density estimate is based on the inverse of the 
volume Wvor,i of the "contiguous" Voronoi cell of each 
datapoint, p(xi) oc l/Wvor,i- The "contiguous" Voronoi 
cell of a point is the cell consisting of the agglomerate of 
all K Delaunay tetrahedra containin^point i as one of its 
vertices, whose volume Wvor,i = 'l2j=i ^Deij is the sum 
of the volumes VdcI.j of each of the K Delaunay tetra- 
hedra. Figure 1 (righthand panel) depicts an illustration 
of such a cell. Properly normalizing the mass contained 
in the reconstructed density field, taking into account the 
fact that each Delaunay tetrahedron is invoked in the den- 
sity estimate at H-M locations, we find at each datapoint 
the following density estimate. 



P(x, 



m(l+M)/Wvor,i 



(2) 
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Having computed these density estimates, we subse- 
quently proceed to determine the complete volume- 
covering density field reconstruction through the linear 
interpolation procedure outlined in Eqn. (1). 

3. Analysis of a cosmological N-body simulation 

Cosmological N-body simulations provide an ideal tem- 
plate for illustrating the virtues of our method. They tend 
to contain a large variety of objects, with diverse mor- 
phologies, a large reach of densities, spanning over a vast 
range of scales. They display low density regions, sparsely 
filled with particles, as well as highly dense and compact 
clumps, represented by a large number of particles. Mod- 
erate density regions typically include strongly anisotropic 
structures such as filaments and walls. 

Each of these features have their own individual char- 
acteristics, and often these may only be sufficiently high- 
lighted by some specifically designed analysis tool. Con- 
ventional methods are usually only tuned for uncovering 
one or a few aspects of the full array of properties. Instead 
of artificial tailor-made methods, which may be insensitive 
to unsuspected but intrinsically important structural ele- 
ments, our Delaunay method is uniquely defined and fully 
self-adaptive. Its outstanding performance is clearly illus- 
trated by Figure 2. Here we have analyzed an N-body 
simulation of structure formation in a standard CDM 
scenario (1^0 = 1, i?o = 50 km/s/Mpc). It shows the re- 
sulting distribution of 128'^ particles in a cubic simula- 
tion volume of 100/i~^ Mpc, at a cosmic epoch at which 
cr{RTH = S/i-^Mpc) = 1. The figure depicts a lO/i-^Mpc 
slice through the center of the box. The lefthand column 
shows the particle distribution in a sequence of frames 
at increasingly fine resolution. Specifically we zoomed in 
on the richest cluster in the region. The righthand col- 
umn shows the corresponding density field reconstruction 
on the basis of the grid-based Triangular-Shaped Clouds 
(TSC) method, here evaluated on a 518^ grid. For the TSC 
method, one of the most frequently applied algorithms, we 
refer to the description in Hockney and Eastwood (1981). 
A comparison with other, more elaborate methods which 
have been developed to deal with the various aspects that 
we mentioned, of which Adaptive Grid methods and SPH 
based methods have already acquired some standing, will 



be presented in Schaap and van de Weygaert (2000) 



A comparison of the lefthand and righthand columns 
with the central column, i.e. the Delaunay estimated den- 
sity fields, reveals the striking improvement rendered by 
our new procedure. Going down from the top to the bot- 
tom in the central column, we observe seemingly compa- 
rable levels of resolved detail. The self-adaptive skills of 
the Delaunay reconstruction evidently succeed in outlin- 
ing the full hierarchy of structure present in the particle 
distribution, at every spatial scale represented in the sim- 
ulation. The contrast with the achievements of the fixed 
grid TSC method in the righthand column is striking, in 



particular when focus tunes in on the finer structures. The 
central cluster appears to be a mere featureless blob! In ad- 
dition, low density regions are rendered as slowly varying 
regions at moderately low values. This realistic conduct 
should be set ofi^ against the erratic behaviour of the TSC 
reconstructions, plagued by annoying shot-noise effects. 

Figure 2 also bears witness to another virtue of the 
Delaunay technique. It evidently succeeds in reproduc- 
ing sharp, edgy and clumpy filamentary and wall-like fea- 
tures. Automatically it resolves the fine details of their 
anisotropic geometry, seemlessly coupling sharp contrasts 
along one or two compact directions with the mildy vary- 
ing density values along the extended direction(s). More- 
over, it also manages to deal succesfully with the substruc- 
tures residing within these structures. The well-known 
poor operation of e.g. the TSC method is clearly borne 
out by the central righthand frame. Its fixed and infiexi- 
ble "filtering" characteristics tend to blur the finer aspects 
of such anisotropic structures. Such methods are there- 
fore unsuited for an objective and unbiased scrutiny of the 
foamlike geometry which so pre-eminently figures in both 
the observed galaxy distribution as well as in the matter 
distribution in most viable models of structure formation. 

Not only qualitatively, but also quantitatively our 
method turns out to compare favourably with respect to 
conventional methods. We are in the process of carefully 
scrutinizing our method by means of an array of quanti- 
tative tests. A full discussion will be presented in Schaap 
and van de Weygaert ( ^000 ). Here we mention the fact 
that the method recovers the density distribution function 
over many orders of magnitude. The grid-based methods, 
on the other hand, only managed to approach the appro- 
priate distribution in an asymptotic fashion and yielded 
reliable estimates of the distribution function over a mere 
restricted range of density values. Very importantly, on 
the basis of the continuous density field reconstruction 
of our Delaunay method, we obtained an estimate of the 
density autocorrelation function that closely adheres to 
the (discrete) two-point correlation function directly de- 
termined from the point distribution. Further assessments 
on the basis of well-known measures like the KuUback- 



Leibler divergence (KuUback and Lciblcr, 1951), an ob- 
jective statistic for quantifying the difference between two 
continuous fields, will also be presented in Schaap and van 
de Weygaert ( |2000| ). Finally, we may also note that in ad- 
dition to its statistical accomplishments, we should also 
consider the computational requirements of the various 
methods. Given a particle distribution, the basic action of 
computing the corresponding Delaunay tessellation, itself 



an 0{N) routine (van de Weygaert, 1991), the subsequent 
interpolation steps, at any desired resolution, are consid- 
erably less CPU intensive than the TSC method (both 
also 0{N)). In the case of figure 2 the Delaunay method 
is about a factor of 10 faster. In the present implemen- 
tation, the bottleneck is Delaunay's substantial memory 
requirement (« 10 x the TSC operation), but a more ef- 
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Fig. 2. A 9-franie mosaic 
comparing the performance 
of the Delaunay density 
estimating technique with 
a conventional grid-based 
TSC method in analyzing 
a cosmological N-body sim- 
ulation. Left column: the 
particle distribution in a 
10/i~^Mpc wide central slice 
through the simulation box. 
Central column: the corre- 
sponding Delaunay density 
field reconstruction. Right 
column: the TSC rendered 
density field reconstruction. 
The colour scale of the den- 
sity fields is logarithmic, 
ranging from 5p/p = — 
2400. 



ficient algorithm will be available in short order. These 
issues will be treated extensively in our upcoming publi- 
cation. 
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The preceding is ample testimony of the promise of 
tessellation methods for the aim of continuous field recon- 
struction. The presented method, following up on earlier 
work by Bernardeau & van de Weygaert ( |1996D , may be 
seen as a first step towards yet more advanced tessellation 
methods. One suggested improvement will be a second- 
order method rendering a continuously differentiable field 
reconstruction, which would dispose of the rather conspic- 
uous triangular patches that form an inherent property of 
the linear procedure with discontinous gradients. In par- 
ticular, we may refer to similar attempts to deal with re- 
lated problems, along the lines of natural neighbour inter- 



polation (3ibson, 1981), such as implemented in the field of 
geophysics (Bambridge et al., 1995; Braun and Sambridge, 
199 



I : • in engineering mechanics (3ukumar, 199S). As 
multidimensional discrete data sets are a major source of 
astrophysical information, we wish to promote such tessel- 
lation methods as a natural instrument for astronomical 
data analysis. 
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